# Code for Table H2: Instrument Validation

rm(list = ls())

## ---------------------------------------
## Load Packages 
## ---------------------------------------


library('AER')
library('ivpack')
library('data.table')


## ---------------------------------------
## Load Data and Functions
## ---------------------------------------
load("../data0812.RData")

directory <- "../functions/"
functions <- list.files(directory)  
loadfunctions <- sapply(functions, FUN = function(x)source(paste0(directory, x)))

## ---------------------------------------
# Construct Instrument
## ---------------------------------------
data0812 <- constructIV(data0812)

## ---------------------------------------
#  Select Last Case Before Election
## ---------------------------------------
data0812 <- lastCase(data0812)


## ---------------------------------------
## Two Stage Least Squares Specifications
## ---------------------------------------

time.controls <- "as.factor(court_time1) + as.factor(court_time2) + as.factor(court_dow) + as.factor(court_shift) + as.factor(totOGS2)"
case.controls <-   "as.factor(any_drug_2) +  as.factor(any_violent_2) + as.factor(fire_arms_2) +  as.factor(any_rob_2) + as.factor(any_dui_2) + as.factor(prior_offender_2)"
demo.controls <- "age_2012 + I(age_2012^2) + Female + as.factor(race) + vote2008 + as.factor(noteli08) + regis_before"

outc.1 <- "vote2012"
endo.1 <- "pti"
inst.1 <- "judgeiv"


form.2 <- formula(paste(inst.1, "~", time.controls))
mod1b <- lm(form.2, data = data0812)
data0812$judgeiv_res <- data0812$judgeiv - predict(mod1b, data0812)

form.1 <- formula(paste(outc.1, "~", demo.controls, "+", case.controls))
m1a1 <- lm(form.1, data = data0812)
data0812$pred_turn <- predict(m1a1)

mod2 <- lm(pred_turn ~ poly(judgeiv_res, 1), data = data0812)

judgeiv_res <- seq(min(data0812$judgeiv_res), max(data0812$judgeiv_res), by = 0.001)
data.pred2 <- data.table(judgeiv_res)

myPredict <- predict(mod2, newdata = data.pred2, interval = "confidence")

pdf("../output/Figures/Figure_H2.pdf", h=8,w=8)  
par(mar = c(6, 6, 4, 6), mgp=c(2,1,0))

hist(data0812$judgeiv_res, 
     freq = F, breaks = 10, las = 1,
     xlab = "Leave-out-case Residualized Pretrial Detention Rate",
     main = "", ylim = c(0, 16))
x <- judgeiv_res
ix <- sort(x,index.return=T)$ix
text(-0.05, 15, paste0("Correlation = ", round(cor(data0812$judgeiv_res, data0812$pred_turn), 3)))

par(new = T)
plot(x, myPredict[, 1], axes=F, xlab=NA, ylab=NA, type = "l", col = "red", ylim = c(0.24, 0.32))
polygon(c(rev(x[ix]), x[ix]), c(rev(myPredict[ ix,3]), myPredict[ ix,2]), col = rgb(0.7,0.7,0.7,0.4) , border = NA)
lines(x[ix], myPredict[ix , 1], col=2, lwd=2)
axis(4, at = seq(from = 0.20, to = 0.40, by = 0.01), 
     labels = seq(from = 0.20, to = 0.40, by = 0.01), 
     col = 'red', col.axis = 'red', las = 1)
mtext("Predicted Turnout Rates", side = 4, col = "red", line = 4)
dev.off()

cat("\nFigure H2: Instrument Validation\n")
cat("Saved Figure H2 in /output/Figures/Figure_H2.pdf\n\n")
